package scu.maqiang.electromagnetics;

import java.util.Arrays;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2T31;
import scu.maqiang.mesh.*;
import scu.maqiang.numeric.Complex;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.IterSCSolver;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SCMatrix;
import scu.maqiang.numeric.SRMatrix;

/**
 * Helmholtz方程求解, 复系数问题
 * @author 马强
 *
 */
public class Muwave {

	public static void main(String[] args) {
		String dir = "muwave\\";
		String fileName = dir + "muwavemesh.msh";
		Mesh2T3 mesh = new Mesh2T3();
		mesh.readMesh(MeshFormat.msh, fileName);
		mesh.toTecplot(dir + "muwavemesh.dat");
		
		int domain1 = mesh.getDomainLabel(new double[]{5.0, 5.0});
		int domain2 = mesh.getDomainLabel(new double[]{17.0, 10.0});
		System.out.println("Two domains: " + domain1 + "\t" + domain2);
		ScalarCFunc coefC = (xy, label, param) -> label == domain2? new Complex(2.0/1.25, 1.0/1.25):new Complex(1.0/1.25, 0.5/1.25);
		FES2T31 fs = new FES2T31(mesh);
		SRMatrix RA = new SRMatrix(fs.GetNdof());
		SCMatrix RC = new SCMatrix(fs.GetNdof());
		fs.assembleStiff(new double[]{1.0}, BVPType.COMMON, RA);
		fs.assembleMass(new ScalarCFunc[]{coefC}, null, BVPType.LUMP_MASS, RC);
		
		SCMatrix LHS = RC.add(-1.0, RA);
//		System.out.println(RA.getElement(609, 609));
//		System.out.println(RA.getElement(609, 605));
//		System.out.println(RA.getElement(609, 546));
//		System.out.println(RC.getElement(0, 0));
//		System.out.println(RC.getElement(0, 1));
//		System.out.println(RC.getElement(1, 0));
//		
//		System.out.println(LHS.getElement(609, 609));
//		System.out.println(LHS.getElement(609, 605));
//		System.out.println(LHS.getElement(609, 546));

		//System.out.println(RA.display());
		Complex[] b = Complex.zeros(fs.GetNdof());
		ScalarCFunc bcFunc = (xy, label, param) -> new Complex(Math.sin(Math.PI * (xy[1] - 15.0)/7.0), 0.0);
		fs.applyBC_MBN(LHS, Direct.All, 1, 2);
		fs.applyBC_MBN(b, Direct.All, new Complex(0.0, 0.0), 1);
		fs.applyBC_MBN(b, Direct.All, new ScalarCFunc[] {bcFunc}, null, 2);
		
		IterSCSolver solver = new IterSCSolver(LHS);

		Complex[] x = Complex.zeros(fs.GetNdof());
		solver.PCGSSSOR(b, x, 1.5, 1);
		//for(int i = 0; i < fs.GetNdof(); i++) {
		//	System.out.println(x[i].image);
		//}
		mesh.toTecplot(dir + "muWaveSol.dat", x);
		
		double[] ff = new double[fs.GetNdof()];
		Arrays.setAll(ff, i -> 1.0e5 * x[i].norm2());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleSource(new double[][] {ff}, BVPType.COMMON, RHS, domain2);
		//System.out.println(Arrays.toString(ff));
		fs.applyBC_MBN(RA, Direct.All, 1, 2);
		fs.applyBC_MBN(RHS, Direct.All, 0.0, 1, 2);
		NewIterSSolver solver2 = new NewIterSSolver(RA);
		double[] xx = new double[fs.GetNdof()];
		solver2.PCGSSOR(RHS, xx, 1.5, 1);
		mesh.toTecplot(dir + "temperature.dat", xx);
		/*
		System.out.println(RC.display());
		System.out.println(RC.getNNZ());
		System.out.println(RC.getElement(0, 0));
		System.out.println(RC.getElement(0, 1));
		System.out.println(RC.getElement(1, 0));
		*/		
	}

}
